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Abstract. Properties of strongly interacting, two-component finite Fermi systems are discussed within 
the recently developed coordinate-space Hartree-Fock-Bogoliubov (HFB) code hfb-AX. Two illustrative 
examples are presented: (i) weakly bound deformed Mg isotopes, and (ii) spin-polarized atomic condensates 
in a strongly deformed harmonic trap. 

PACS. 21.60.Jz Nuclear Density Functional Theory - 31.15.Es Atomic Density Functional Theory 
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1 Introduction 

Superconductivity and Cooper pair formation are generic 
features of strongly-interacting many-body Fermi systems. 
In the context of the Density Functional Theory (DPT), 
the Hartree-Fock-Bogoliubov (HFB or Bogoliubov-de-Gen- 
nes) framework has been widely used to treat pairing cor- 
relations in nuclei (see, e.g., [TJOE]) and ultracold atom 
gases (see, e.g., 3',T,'5','6',T). The superiority of the HFB 
method over the conventional BCS approximation becomes 
particularly apparent in the context of weakly bound sys- 
tems, such as drip-line nuclei, where the coupling to the 
scattering continuum becomes essential T . 

The HFB equations can be solved in several ways (see 
Ref. |8j for a recent overview) . In the commonly used con- 
figuration space approach, the quasi-particle orbitals are 
expanded in a suitable single-particle basis. A number of 
HFB codes were developed by employing the harmonic 
oscillator (HO) eigenstates [2l|9]. However, use of the HO 
basis is questionable in the limit of both weak binding and 
very large deformations, which require the use of unreal- 
istically large configuration spaces to guarantee conver- 
gence. In both situations, the coordinate-space approach 
to the HFB problem p1llOilll| is superior. 

A number of coordinate-space techniques have been 
developed over the years, and their performance strongly 
depends on the size and self-consistent symmetries of the 
spatial mesh employed [T2l[T3]. While the direct itera- 
tive diagonalization procedure in the coordinate space is 
computationally intensive, the advent of teraflop super- 
computers enables us to carry out large-scale DFT com- 
putations of complex physical systems in non-spherical 
spatial boxes. The recently developed parallel 2D-HFB 



solvers utilizing the B-spline technique offer excellent ac- 
curacy when describing deformed weakly bound nuclei |14[ 
[8] . Solving the HFB equations in a 3D coordinate space is 
not a simple task; it is worth noting that several develop- 
ments are underway, such as a general-purpose 3D-HFB 
solver based on multi-resolution analysis and wavelet ex- 
pansion [T?[|E]. 

We recently released a 2D coordinate-space code HFB- 
AX that solves the HFB problem using B-splines [8j. Its 
high precision has been explicitly demonstrated by testing 
against the HO basis expansion method, wavelet method, 
and other HFB codes. In this work, we apply hfb-AX to 
two problems of current interest. First, we study the drip- 
line Mg isotopes using the SLy4 TS] energy density func- 
tional. The heaviest-known ""^Mg isotope has recently been 
produced [17] , and is expected to be weakly deformed [18] . 
Here, we systematically compare the differences between 
HO expansion and coordinate-space HFB calculations. 

In addition to nuclear calculations, hfb-AX has re- 
cently been applied to cold atomic systems within the 
Superfluid Local Density Approximation (SLDA) [TlfTO]. 
Studies of pairing properties of polarized cold atom gases 
are indeed of considerable experimental [2(1, 21] and theo- 
retical [UlUIS] interest. The separation between paired and 
normal phases has been observed in a gas with unequal 
numbers of two spin components trapped in an extremely 
deformed potential [20 . Calculations indicate that experi- 
mental results depend on the number of fermions and trap 
asymmetry Theoretical simulations based on the HO 
basis method are, however, limited to spherical shapes or 
small deformations [11[5]. In this work, we calculate the 
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Table 1. Comparison between hfb-AX and hfbtho with SLy4 p-h functional and mixed pairing for neutron-rich even-even 
nuclei ^"""''Mg. AU energies are in MeV. The quadrupole moments Q20 are in fm^. See text for more details. 

Nuclei HFB-AX HFB-AX HFBTHO HFB-AX HFBTHO HFB-AX HFBTHO HFB-AX 



Q2O Etot Etot Ej^in Ef^in Epair Epair 



-Mg 


-fl48.0 


-3.63 


-257.49 


-257.48 


553.51 


553.54 


-5.28 


-3.98 


^'^Mg 


+175.4 


-2.49 


-263.02 


-263.06 


585.89 


585.63 


-0.02 


-0.12 


3«Mg 


-hl87.9 


-1.37 


-265.40 


-265.46 


626.39 


625.16 


-2.49 


-2.59 


*°Mg 


-F221.6 


-0.45 


-266.96 


-267.13 


655.66 


653.98 








''Mg 


-136.1 


-0.22 


-264.83 


-264.91 


699.68 


695.00 


-5.11 


-4.34 


''Me 


-107.8 


+0.15 


-264.19 


-264.56 


730.35 


723.02 


-2.68 


-3.44 



properties of a Fermi gas of 1000 polarized atoms trapped 
in an extremely elongated HO potential. 



2 Nuclear calculations 

In the coordinate-space representation, the HFB equations 
can be written as [HE]: 



""^y Mra,rV) 



^(2)(i.'ct') 



h{r(7, r'a') 
h{ra, r'a') + A 



(1) 



E 



where (r, a) are the particle spatial and spin coordinates, 
h{ra,r'a') and h{ra,r'a-') are the particle-hole (p-h) and 
particle-particle (p-p) components of the HFB Hamilto- 

nian, respectively, ipn\ra) and Tpn\ra) are the upper 
and lower components of the HFB wave function, and A 
is the chemical potential. The spectrum of quasiparticle 
energies E is discrete for \E\ < —A and continuous for 
\E\ > —A. By requiring that the eigenfunctions vanish at 
the edge of the box (box boundary conditions) , the parti- 
cle continuum becomes discretized. 

In the axial geometry, the third component of the single- 
particle angular momentum, i7, is a good quantum num- 
ber. The HFB wave function can thus be written as ^^^{r) 
where r — z), q—±^ denotes the cylindrical isospin 
coordinates, and i7=±i, ±|, ±|, . . .. We also assume that 
the reflection symmetry is conserved; hence the parity tt is 
a good quantum number. Consequently, the Hamiltonian 
matrix becomes block diagonal, which enables paralleliza- 
tion. 

In nuclear calculations, the p-h channel is often mod- 
eled with the Skyrme energy density functional, while a 
zero-range isovector pairing interaction is often used in the 
p-p channel. By using the zero-range interactions, Eq. ([1]) 
becomes local and easier to solve. In this study, we used 
the SLy4 [16j Skyrme functional and the mixed 5 pairing 
[22] . The pairing strength has been fitted to reproduce the 
average neutron pairing gap in ^^"^Sn [9]. [The numerical 
effort depends on the box size (pmax, ^max), the largest dis- 
tance between neighboring mesh points in the grid h (the 
B-spline grid is not uniform), and the order of B-splines 



M]. We demonstrated [8] that /i=0.6fm and M—Xi guar- 
antee excellent precision of calculations. For heavy elon- 
gated nuclei, such as those on the way to fission, large 
2D boxes are required. To accelerate convergence of self- 
consistent iterations, we employed the modified Broyden 
mixing [23] which is significantly faster than the standard 
linear mixing method. In the present calculations for the 
Mg isotopes, we used a box of Pmax=Zmax=^^ fni. 

Table I displays the results of calculations for •^'^^'''^Mg 
isotopes with hfbtho [2] and HFB-AX. Our hfbtho cal- 
culations were carried out in a configuration space of 20 
HO shells. We compare the total binding energy Etoti ki- 
netic energy Ekim pairing energy Epair, total quadrupole 
moment Q2o, and neutron chemical potential A„. We can 
see that while the two codes yield very similar total bind- 
ing energies, their differences increase towards the neu- 
tron drip line. For example, in ^"^Mg the difference in Etot 
is only lOkeV between the two codes, while it becomes 
170 keV in ^°Mg, As discussed in [5], the kinetic and pair- 
ing energies are slightly different in hfbtho and hfb- 
AX because of different continuum discretization: hfbtho 
predicts kinetic energies that are systematically larger than 
in HFB-AX, especially for nuclei near drip lines. 

In our calculations with SLy4, "'"Mg is the last Mg iso- 
tope that is stable against the two-neutron emission. This 
is consistent with recent experimental evidence [17]. Ac- 
cording to Table I, its ground state is predicted to have a 
well-deformed, prolate shape. (A deformed halo structure 
in "^^Mg has also been predicted in Rcf. [24J.) Interestingly, 
due to deformed shell gaps, static proton and neutron pair- 
ing vanishes in this nucleus. The heavier even-even isotope 
of ^^Mg is predicted to have a negative neutron chemi- 
cal potential but is unbound to two-neutron emission. As 
this nucleus is expected to have appreciable oblate defor- 
mation, however, the particle stability of ^^Mg could be 
enhanced against the two- neutron decay [18]. 



3 Atomic calculations 



The strongly interacting, polarized Fermi gas can be de- 
scribed by introducing mismatched chemical potentials, 
A| and A|, corresponding to spin-up (majority) and spin- 
down (minority) states, respectively. The HFB equations 
of the resulting Two-Fermi Level Approach (2FLA) read 
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A*{r) -Ho{r) + xJ U,(r) ) ' ^nv^{r) 



As usual, we define the average chemical potential A — 
(Af + A|)/2 and the difference 2As = A| — A|. The diago- 
nalization of Eq. ^ is equivalent to Eq. ([T|) by replacing 
the two chemical potentials with A, and adding a cranking 
term involving spin projection. As discussed in Ref. [3j, 
the 2FLA is equivalent to the standard blocking proce- 
dure, and Xg can be viewed as a rotational frequency that 
generates spin polarization. 

The associated polarization density m(r) — p^(r) — 
Pl{r), total density p{r) — pi{r) + Pi{r), pair density 
K(r), and pairing potential A{r) are: 



m{r) 


- E 


(|u,(r)p + K(r)p), 


(3) 




0<Ei<\s 






pir) 


= m{r) + 




(4) 










/v(r) 






(5) 










A{r) 


= -9eff{r 


)K(r), 


(6) 



where geff(r) is the regularized pairing strength. The reg- 
ularization procedure is introduced because the kinetic 
and pairing energies diverge with energy for zero-range 
pairing forces 19, 2 5^. 




The energy density functional of SLDA (in units fi^m^l) 
is given by [7j: 



£{r) 



r(r) 



■/3 



3(37r^)^/3p5/3(^) 
10 



9eff 



pi/3{r) 



(7) 



where a and (3 are dimensionless parameters. In the present 
calculation, we took the axial external HO trap V{p, z) = 
^muP'{p^ + /rf'). The parameter r] defines the trap's an- 
isotropy and we took lo=\. In the experiment of Ref. |20j . 
an extremely large deformation of r/ ~50 was employed. 
Previous theoretical work based on HO expansion con- 
sidered either a spherical geometry [5] or a fairly small 
value of rj <4 j4]. By taking advantage of the coordinate 
space formalism, we could reach very large deformations of 
?7=10. Unlike in nuclei, the spin-orbit coupling is absent 
in the atomic problem. Therefore, the dimension of the 
Hamiltonian matrix for the gas is half of the nuclear one, 
and very large spaces can be reached. We adopted a rect- 
angular box with dimensionless Pmax—W and z,„aa;=50 for 
the calculation of 1000 atoms (A^=iV|+A^x=1000). The 
polarization, P—{N^ — Ni)/N, is realized by adjusting 
the value of A^. The resulting quasi-particle atomic spec- 
trum is much denser than a typical nuclear one. A finer 
mesh size h=OA than in nuclear calculations, and M=13 
order B-splines guarantee convergence of our SLDA cal- 
culations. 
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Fig. 1. (Color online) The calculated density contours of Fermi 
gas with a spin polarization of P=0.21: majority spin state 
density Pi{r) (a); minority spin state density Pi{r) (b); polar- 
ization density (c); and the pair density K(r) (d). The panels 
are plotted with dimensionless length of 80 and width of 20. 
The corresponding density values are given in the legend. To 
facilitate comparison, k is multiplied by a factor —0.5. 



Fig. 2. (Color online) The polarization (a) and pair (b) densi- 
ties of trapped fermions at different polarizations P. The main 
effect of the polarization is to decrease the superfluid core size. 

Figure [T] shows the density contours for the case of 
214 unpaired particles (P=0.21). The distribution of pair 
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density K{r) coincides with that of Pi{r). That is, the 
spin-polarized (unpaired) fermions shown in Fig. [ijc) are 
distributed around the paired superfluid core. This spa- 
tial separation between normal and superfluid phases is 
similar to what has been seen experimentally (20j . 

The densities m^p = 0, z) and k{p = 0, z) correspond- 
ing to different polarizations P are displayed in Fig.[2l The 
effect of phase separation is clearly seen. As discussed in 
Ref. [5] , the main effect of polarization is to decrease the 
size of the superfluid core. The superfluid-to-normal tran- 
sition is not sharp, as expected in the flnite system. In the 
intermediate region of small (but nonzero) a and nonzero 
m, the pair density exhibits small-amplitude oscillations. 
As discussed in Refs. [4j[5], such behavior is characteris- 
tic of Fulde-Ferrel-Larkin-Ovchinnikov (FFLO) phase [26] 
associated with the magnetized superfluid. 

4 Summary 

We used the recently developed coordinate-space HFB 
code HFB- AX to study the strongly interacting superfluid 
Fermi systems, including nuclei and cold atomic gases. 
We first calculated the properties of weakly bound de- 
formed Mg isotopes by solving Skyrme HFB equations. 
We conclude that the two-neutron drip line in the Mg 
chain predicted with SLy4 density functional corresponds 
to ^°Mg, but ''^Mg can be long-lived due to shape coexis- 
tence effects. Secondly, we investigated the density distri- 
butions of a polarized cold atomic gas in an extremely 
deformed external trap. The calculated density profiles 
clearly show the separation between the paired and po- 
larized normal phases, in agreement with experimental 
findings. In the intermediate magnetized superfluid re- 
gion, the pairing field exhibits small amplitude oscilla- 
tions, consistent with FFLO behavior. We conclude that 
the coordinate-space HFB framework is a very useful tool 
for the description of weakly bound and/or extremely de- 
formed and polarized superfluid Fermi systems. 
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DE-FC02-07ER41457 with UNEDF SciDAC Collabora- 
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tional Center for Computational Sciences at Oak Ridge 
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